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ABSTRACT 

Close-in planetary systems detected by the Kepler mission present an excess of periods ratio that are just slightly larger than some 
low order resonant values. This feature occurs naturally when resonant couples undergo dissipation that damps the eccentricities. 
However, the resonant angles appear to librate at the end of the migration process, which is often believed to be an evidence that the 
systems remain in resonance. 

Here we provide an analytical model for the dissipation in resonant planetary systems valid for low eccentricities. We confirm that 
dissipation accounts for an excess of pairs that lie just aside from the nominal periods ratios, as observed by the Kepler mission. In 
addition, by a global analysis of the phase space of the problem, we demonstrate that these final pairs are non-resonant. Indeed, the 
separatrices that exist in the resonant systems disappear with the dissipation, and remains only a circulation of the orbits around a 
single elliptical fixed point. Furthermore, the apparent libration of the resonant angles can be explained using the classical secular 
averaging method. We show that this artifact is only due to the severe damping of the amplitudes of the eigenmodes in the secular 
motion. 
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1. Introduction 

Dissipation due to tidal interactions is a possible mechanism 
that explains the abundance of planetary systems that lie near 
but not at exact mean-motion commensurability. Papaloi zou &| 
|Terqu em ( 20101 (in the case of three planets Laplace resonances) 
and |Papaloizou| ( [201 1\ (for two planets resonances) showed that 
planets that have been temporarily locked in resonance due to 
differential migration could have their periods ratios to depart 
from strict commensurability due to the circularization of their 
orbits by tidal interactions with the star. More recently, Batygin 
|& Morbidelli1 ( |20T2l i and |Lithwick & Wu| ( |20T2l i use similar ef- 
fects to explain the excess of systems of two planets that lie near 
resonances but with planets slightly farther from each other than 
the nominal mean-motion commensurability ratio. 

One of the most intriguing features that is common in these 
different studies is the observation that resonant angles continue 
to librate far from exact commensurability and the authors leave 
unanswered the question of determining if these systems are in 
resonance or not. It seems important to clarify this point and 
to understand why resonant angles can librate so far from exact 
commensurability. 

Our analysis is based on the study of the phase space of two 
planets in mean-motion resonance (MMR) in the conservative 
case (without any dissipation). This is the object of section|2] We 
pay a particular attention to apsidal corotation resonances (ACR) 
which play a major role in the dynamics of these systems and in 
the understanding of the topology of the phase space. ACR have 
been extensively studied both in the asteroidal restricted problem 
(e.g. Ferraz-Mello et al.| 1993 1 and the planetary problem (e.g. 
Hadj i demetriou|2002||Michtchenko et al.|2006] > . 



Most of the studies on the subject have been made using nu- 
merical (or semi-analytical) models that remain valid for arbi- 
trary values of the eccentricities but which do not always pro- 
vide a global picture of the dynamics. In the present study we 
are only concerned in the dynamics at low eccentricities since 
our aim is to understand the motion at the end of the circulariza- 
tion process. A completely analytical model is thus well suited 
in this case. Analytical studies of planetary MMR have already 
be done up to degree two in eccentricities in the cases of the 2: 1 
( [Callegari et al.|2004) > and the 3:2 ( |Callegari et al.|2006| > MMR. 

The dissipative case (studied in section [3]) is modeled us- 
ing the conservative case as a basis and very simple and general 
prescriptions for the dissipation. Our study is mainly aimed at 
understanding the impact of tides on the dynamics of resonant 
planets pairs and, in this case, we follow the prescriptions intro- 
duced by Papaloizou] ( |201 1[ ). However we show that the differ- 
ential migration process that allows a resonant locking of both 
planets can also be accounted for in our model. This process 
has already been widely studied (e.g.|Lee & Peale 2002 ; F erraz-| 
Mello et al.|[2003| |j~xe1|2"004"t |Beauge et al.l|2006), and is also 



considered in Papaloizou & Terquem ( 2010 ); Papaloizou ( 201 1 



Eventually, we treat this perturbation of the conservative case by 
following the lines of Laskar et al. (2012 1. 



In section [4] we show that the final state of the resonant sys- 
tems that undergo a circularization process is very well charac- 
terized by the secular normal form. We explain, with the secular 
problem, why resonant angles appear to librate even far from 
resonances. 

Finally, we present in section|5]the results of two numerical 
simulations that confirm and illustrate the different mechanisms 
that we highlight with our analytical model. 
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2. Dynamics of two resonant planets in the 
conservative case 

2.1. Model 

We study in this section the case of two planets orbiting a star in 
the same plane and without any dissipative force. This problem 
is usually referred to as the planar planetary three body problem. 
We note with a subscript 1 the internal planet and 2 the external 
one. The star is referred to as body 0. Masses are noted w,. For 
both planets we define: fij = @(mo+nii), and/?,- = mom,7(mo+m;), 
where @ is the gravitational constant. We note a, the semi-major 
axis of planet i, and e, its eccentricity. 

We use Delaunay canonical pairs of astrocentric coordinates 
for both planets. The actions are the circular angular momentum 
and the angular momentum of both planets: 



A 



(1) 

(2) 



The associated angles are the mean anomalies (M,) and the ar- 
guments of periastron («,■) of both planets. 

The Hamiltonian of the system reads, in these coordinates: 



K = 9* (A) + <H X (A, G,M,cS), 



(3) 



where 'Hq is the Keplerian part and r H[ is the perturbative part 
due to planet-planet interactions taking into account both direct 
and indirect effects. The Keplerian part depends only on A, (i = 

1,2): 



Ho 



2 ,,2 

-z 



1=1 



2A 2 ' 



(4) 



Whereas the perturbative part depends on all eight Delaunay co- 
ordinates. We do not need to express the explicit form of 'Hi at 
this point but it could be seen as a Fourier series of all four angles 
(with coefficients depending on actions). 

At first sight, we have to deal with a four degrees of free- 
dom differential system. However, we will now introduce two 
well-known transformations which will reduce this problem to 
two degrees of freedom. The first reduction makes use of the 
total angular momentum conservation and is always valid. The 
second one corresponds to an averaging of the equations of mo- 
tion and is only valid near a specified mean-motion resonance 
(MMR). More precisely, it is only valid far from all other MMR. 

2.1 .1 . Conservation of the total angular momentum 

In the Delaunay coordinates system, the total angular momen- 
tum conservation (G = G\ + G 2 ) is equivalent to the fact that the 
Hamiltonian depends only on the value of the difference of argu- 
ments of periastron: Aco = cj] - a> 2 and not on individual values 
of both angles. The reduction resulting from the conservation of 
the angular momentum is performed by using the coordinates 
Au>, 0J2 instead of u\, co 2 . These two new angles are respectively 
canonically conjugated to the actions G\, G. Therefore, the sys- 
tem of coordinates (M\ ,A\, M 2 , A 2 , Aoj, G\,a)2,G) is canonical 
and with these coordinates, C02 does not appear anymore in the 
expression of the Hamiltonian. G is a first integral of the equa- 
tions of motion. We can thus study the three degrees of freedom 
system constituted of (Mi, Ai, M 2 , A 2 , Aco, G\) and consider G 
as a parameter of this system. 



2.1 .2. Averaging the Hamiltonian near a MMR 

We now suppose that the system is near a mean-motion com- 
mensurability of the form (p + q):p, that is: 

- pM x + (p + q)M 2 ~ 0. (5) 

We may then introduce the argument of the (p + q):p MMR: 



P 



M 2 , 



(6) 



o~ - --Mi 
1 

such as & « . 

We then construct the resonant normal form to the first order 
of the planets masses (with respect to the stellar mass) by av- 
eraging over all rapid angles, i.e. all combinations of the mean 
anomalies that are not harmonics of cr. By doing this averaging 
we need to introduce a new set of variables which corresponds 
to the averaged problem. The change of coordinates is close to 
identity so we can identify both sets of coordinates in first ap- 
proximation. Strictly speaking the transformation to the initial 
set of coordinates reintroduces high frequencies (e.g. Mo rbidelli| 
2002). Before constructing the normal form, we introduce the 
change of variables: 



q 



a, 



(7) 



(8) 



F = ll + -^A 1+ A 2 , 

and we consider the canonical set of coordinates: 
(cr, I, M2,r, Aco, Gi). With this set of coordinates, the reso- 
nant normal form is obtained by averaging over M%. Thus, 
by definition, M2 does not appear anymore in the averaged 
Hamiltonian, and T is a first integral of the averaged motion. 
Therefore, the system is reduced to a two degrees of freedom 
problem with four canonical coordinates : (cr, /, Aw, G 1 ) and two 
parameters: G and T. However this set of coordinates is not 
very well suited for the study of low-eccentricities systems and 
it is a lot more convenient to introduce rectangular canonically 
conjugated variables similar to Poincare variables. 

2.1.3. Rectangular coordinates 

We first introduce the two resonant angles which correspond to 
both planets (e.g.|Ferraz-Mell o et al.|1993] l: 



a" i 



-(-5) 



Aa> = --A x + (l + -]A 2 -co\, 



(9) 



cr 2 = lt--Aco = --Ai + \\ + -\a 2 -co2, (10) 
q q \ q) 

where /I, = M, + w, is the mean longitude of planet 2. We com- 
plete this change of variables with the canonically conjugated 
actions: 

f P : 



1 



-l-Gi = Ai - Gi, 



(11) 



(12) 



/ 2 = 1 1 + C y + Gj = A 2 - G 2 + G - T. 

Note that, since G and F are constants of the motion, we can 
redefine the action I2 as: 

/ 2 =A 2 -G2, (13) 

without any change for /1, cr] and <r 2 . 

The associated rectangular coordinates are then: 



(14) 
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2.1 .4. Elimination of the first integral r 

The dynamics of our two degrees of freedom system should de- 
pend on the value of both first integrals: F and G. However, it is 
possible to eliminate the dependency in T by dividing all actions 



by its value: 




A A '- 

Ai = T • 


(15) 


c - di 
- y> 


(16) 


G 

G = - =Gi+G 2 , 


(17) 


h = ^=Ai-Gi, 


(18) 


h = y - A2 _ G2, 


(19) 


Xi 

Xi = — . 

Vr 


(20) 



With these renormalized variables, Hamilton equations are still 
valid if the Hamiltonian is also divided by T: 



r ' 



(21) 



However, the Hamiltonian < 7if, and thus the dynamics, still de- 
pends on the value of F. The complete elimination of T is 
achieved by renormalizing both energy and time scales: 



<H = F <H = ra, 
t 

T= T3- 



(22) 
(23) 



It can be verified that 'H does not depend anymore on F and 
Hamilton's equations now reads: 



dxj 
dr 



.an 

dxi 



(24) 



This means that the value of F does not influence the dynamics 
of the system except by changing the scales of distance, time 
and energy (respectively by a factor of F 2 , F 3 and 1/F 2 ). With 
this renormalization, we reduced the number of parameters of 
our system from two (G, F) to one (G = G/T). 

2.1 .5. Explicit form of the Hamiltonian 

We need to express the Hamiltonian *H as a function of x,- 
(i = 1,2) and G. Let us start with the Keplerian part. The 
renormalized Keplerian part CKo) has the same form as "Ho (see 
Eq. Q) if we replace A, with A,-. Moreover, A, are simple func- 
tions of Xj and G: 



A, 



(G + D-l), 



A 2 = 1-ll + i )a 1 = [i + P -)(G+D)- P 



(25) 
(26) 



where T> is the (renormalized) angular momentum deficit 
(AMD, e.g. |Laskar|2000) defined by: 



D = Y J x i x i (=A 1+ A 2 -G). 



(27) 



Therefore, the Keplerian part is obtained by substituting A, 
values in the expression of fio. This expression is then expanded 
to a given degree in eccentricities (i.e. to a given power in Xi vari- 
ables). Note that the Keplerian part is now expressed as a power 
series of eccentricities even if it only depends on semi-major 
axes in Delaunay coordinates system. It should be remarked that 
Xi variables only appear in the Keplerian part with T> (which is 
of degree two and symmetrical for both planets). As a conse- 
quence the Keplerian part only contains terms of even degree in 
eccentricities. 

For the perturbative part, we use the expansion method pre- 
sented in Las kar & Robutel] (|1995|l and impleme nted in the 
algebraic manipulator TRIP dGastineau & Laskar 2011). The 
Poincare coordinates used in [Laskar & Robutel| ( |1995| i7 noted 
x, and x', are related to our rectangular coordinates with the re- 
lations: 



Xi e" 
x 2 e u 



with: 



- P -A 1+ [l + P -W 



(28) 
(29) 



(30) 



Laskar & Robutel ( 1995 ) method is aimed at computing the 
perturbative part of the Hamiltonian in power series of x, x' and 
in Fourier series of the mean longitudes A\, Xi (with coefficients 
depending on A,). The averaging method described in section 
2. 1.2 consists in selecting in this expansion the terms that do not 



depend upon combinations of the mean longitudes other than ctq 
(and its harmonics). 

Then it is straightforward to express 'K; as a function of x; 
and G, by substituting x, x' and A, by their values and renor- 
malizing by T. We obtain a power series of x, with coefficients 
depending on G. 

The ratio between the perturbative and the Keplerian parts is 
of the order of the mass ratio between the planets and the star. To 
be consistent, the Keplerian part must be expanded to a higher 
degree in eccentricities than the perturbative part. 

As an example, for a first order resonance (of the form 
(p + l):p), the Hamiltonian expanded to the fourth degree in 
eccentricities for the Keplerian part and to the second degree for 
the perturbation (noted degree 4-2) has the form: 



<H = K (1) D + K {A) D 2 



S^D+Sfx^ 



? (2) 



*l) ' 



+ s 2 

R?(x 2 



(2) - , c (2), - , - x 



R^\x\ + x\)+Rf{x: 



x 2 ) 

2 ) +R?hxix 2 + x{x 2 ), 



(31) 



!=1 



where all coefficients (K, S , R) depends on G and on the masses, 
and K stands for Keplerian terms, S for secular terms and R for 
resonant terms (see appendix [A] for a description of their com- 
putation). For a resonance of order two (or more), the expres- 
sion would be similar but there would not be resonant terms of 
the degree one. These terms are very important in the dynamics 
since they introduce constant terms (non-zero right-hand side) 
in Hamilton's equations (Eq. (|24|). When there are no first de- 
gree terms (resonance of order two or more), x, = (2 = 1, 2) is 
always a fixed point. On the contrary, for first order resonances, 
this fixed point does not exist and an initially circular system will 
not stay circular. Note that since the Hamiltonian is expanded as 
series of eccentricities, this model is only valid for low eccen- 
tricities. 
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2.2. Study of the dynamics 
2.2.1 . Energy levels 

We consider a first order resonance and construct the resonant 
normal form developed up to degree 4 for the Keplerian part and 
to first degree for the perturbation (noted degree 4-1). For given 
values of the masses of the three bodies, a given mean-motion 
commensurability (p + q): p, and a given value of G, we are able 
to compute the needed coefficients K, S and R. We can then look 
at the energy levels of this Hamiltonian. Since our problem has 
two degrees of freedom and four dimensions, we cannot have a 
global view of these energy levels but we can represent them on 
section planes. Figure [T] shows the energy levels, represented in 
different section planes (see the legend), in the case of the 2: 1 
MMR with a star mass of niQ = m Q (solar mass), planets masses 
of m.\ = m.2 = m m (Earth's mass), and with G = Go(l - 1CT 4 ). Go 
is the value of the total circular angular momentum at nominal 
resonance. At nominal resonance, the exact commensurability of 
mean-motions induces the relation: 



«i 



n 2,r P \M,rj 



mi) Pi 



. 1/3 

mo + m% \ mi 



niQ + mi 



ni2 



(32) 
(33) 



Together with the relation (l + 2)A 1 , r + A 2 , r = l(seeEq. (g), 
this imposes the values Ai r and A 2r at nominal resonance: 




-O.02-O.01 0.01 0.02 

e\ cos (Ti 



-0.02-0.01 0.01 0.02 
ei sin cft 



Fig. 1. 2:1 resonance energy levels in section planes defined 
by: <?2 = (top left), cos cr\ = coscr 2 = (top right), 
sin en = sin <x 2 = (bottom left), and e x — (bottom right). 
Stable ACR solutions are highlighted with colored dots (red 
and green). Unstable ones are highlighted with colored crosses 
(blue). The star mass is given by mo = m B and planets masses by 
ni\ — ni2 = m e . The constant G is set to G = Go(l - 10~ 4 ). The 
Hamiltonian is developed up to degree 4-1. 



Al,r = 
A 2 ,r = 



1 + 1 + 



2/3 



-1 



(34) 
(35) 



Go is then given by: 



G = A v + A 2 , r = { 1 + 



1/3 



(36) 



On Fig.[T] we distinguish three areas in the different section 
planes (except in the plane (sin <ji , sin <x 2 ) on the top right of the 
graph). There are two zones of circulation: internal circulation 
for low eccentricities and external circulation for high eccentric- 
ities. Between these two circulation areas, we observe a libration 
zone (banana-shaped level curves) separated from both circula- 
tion areas by two separatrices. Fixed points of the system (see 
next section) are marked with colored dots (for stable ones) and 
crosses (for unstable ones). The red dot corresponds to the libra- 
tion center in Fig.[T] bottom-left. 

Figure[2]shows the energy levels on the section plane defined 
by sino - ! = sin<x 2 = for increasing values of the parameter G 
and for the 2:1 and the 3:2 MMR. The masses are the same as 
for Fig. JT](mo = m Q , mi - m 2 = m m ). Figure |2]C1 is the same 
as Fig. [lj bottom-left. When G decreases (Figs. |2]A1,B1), the 
libration area moves to higher eccentricities and the internal cir- 
culation area takes more space. On the contrary, for higher values 
of G (Figs.|2jDl to Fl), the internal circulation and the libration 
areas tend to shrink and eventually completely disappear leaving 
only the external circulation area. 

We observe the same evolution for the 3:2 MMR (Figs.|2jA2 
to F2) and it seems that this behavior is common to all first order 
MMR. 



Even if these section planes give an insight on the nature of 
the motion, it is important to keep in mind that the real motion 
occurs on the total four dimensional phase space and that these 
sections cannot provide a global picture of the dynamics. The 
easiest way to have a good insight on the global structure of the 
phase space is to look at the positions and natures of the fixed 
points of the system. These fixed points are commonly referred 
to as apsidal corotation resonances or ACR (e.g. |Ferraz-Mello 
let al.|1993|). 



2.2.2. Fixed points 

The positions of the fixed points are obtained solving Hamilton's 
equations (Eq. (|24|). For instance, at degree 4-2 (Eq. (31 1): 



= (K <2) + 2K^D + S (0) +Sf } ) xi +Sf}x 2 
= S™X! + (# (2) + 2K^D + S (0) +Sf)x 2 



+Rf 2 ] xi +2Rfx 2 +R 2 l) . 



(37) 



(38) 



This imply to solve a system of four polynomial real equations 
(real and imaginary parts of Eqs. p7| , ( |38] >) with four real un- 
knowns (real and imaginary parts of x\, x 2 ). We solve it using 
the Maple RootFinding[Isolate] function (see [Rouillier 1999| l. 
Then, for each solution, we linearize the equations of motion 
around the fixed point in order to compute the eigenvalues and 
to distinguish between elliptic (purely imaginary eigenvalues) 
and hyperbolic (nonzero real parts) fixed points. 



Degree 4-1 

Figure [3] shows the positions and nature of ACR solutions as 
functions of G, in case of the 2:1 MMR with the same masses 
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-0.06-0.04-0.02 0.02 0.04 0.06 
e\ cosci 



-0.06-0.04-0.02 0.02 0.04 0.06 
e\ C0S(Ti 



-0.02 -0.01 0.01 0.02 
ei coscti 



0.02 



0.02 



0.01 



0.01 



0.112 



-0.02 



0.02 



-0.02 




-0.02 -0.01 0.01 0.02 
e\ coso*i 



-0.02 -0.01 0.01 0.02 
e\ cos (j\ 



-0.02 -0.01 0.01 0.02 
e\ cos (T[ 



0.02 



0.02 



0.02 



0.02 



for the 2:1 resonance (Al to Fl), and the 3:2 resonance 



Fig. 2. Energy levels sections in the plane defined by sin en = sincr 2 

(A2 to F2). The constant G/G Q - 1 is set to -1(T 3 (Al, A2), -5 x 1(T 4 (Bl, B2), -1(T 4 (CI, C2), -3.2 x 1(T 3 (Dl, D2), -3 x 10 
(El, E2), and +1CT 3 (Fl, F2). Stable ACR solutions are highlighted with colored dots (red and green). Unstable ones are highlighted 
with colored crosses (blue). The star mass is given by mo - m Q and planets masses by ;«i = mi = »%. The Hamiltonian is developed 
up to degree 4-1. Note that the scales are different for graphs Al, A2, Bl, and B2 than for other graphs. 



than before. The Hamiltonian is still developed up to degree 4- 
1 . On this figure, continuous lines represent stable fixed points 
and dashed lines represent unstable ones. The choice of colors 
is consistent with Figs. [T] and [2] We gave the positions of ACR 
only in the directions e\ coso - ! and e2 cos o"2 because all ACR 
solutions that we found have sincri = sino^ = 0. On the left 
of Fig. [3] (lowest values of G), there are two elliptic and one 
hyperbolic ACR. This corresponds to the situation of Figs.[T]and 



[2]A1 to Dl. The red elliptic ACR corresponds to the center of 
the libration area. The green one corresponds to the center of the 
internal circulation area. The blue hyperbolic ACR corresponds 
to the crossing of internal and external separatrices. 

Around G = Go, we observe a bifurcation and the green and 
the blue ACR disappear whereas the red one remains the only 
one to survive and it stay stable. This corresponds to the situation 
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0.15 




0.1 




0.05 


b 




o 





o 






-VJ.VJJ 




-0.1 




-0.15 




0.04 




0.03 




0.02 




0.01 






O 





o 




cs 


-0.01 




-0.02 
-0.03 
-0.04 

-0.0025 -0.002 -0.0015 -0.001 -0.0005 
G/Go - 1 



0.0005 



Fig. 3. 2:1 ACR positions as functions of G. The star mass is 
given by toq = m e and planets masses by ni\ = m-i = m @ . The 
Hamiltonian is developed up to degree 4-1. Elliptic (stable) ACR 
are plotted using continuous lines whereas hyperbolic (unstable) 
ACR are plotted using dashed lines. At this degree of develop- 
ment, all ACR have sin cr\ = sin cr^ = 0. 




0.0001 



Fig. 4. 2: 1 ACR positions in the plane (G, a) in the same condi- 
tions as for Fig. [3] 



observed in Fig.|2jEl,Fl and which cannot be considered as a 
resonant situation. 

In order to have a better view of the spatial positions of these 
ACR we plotted on Fig. [4] the semi-major axes ratio of ACR 
solutions as functions of G. We observe that the center of li- 
bration corresponds to values of a lower than the nominal res- 
onant value. This is equivalent to values of PiJP\ greater than 
the nominal resonant value (2). It means that planets are farther 
away from each other than the nominal resonance distance. On 
the contrary, the hyperbolic ACR and the center of the internal 
resonant area correspond to higher values of a, thus two planets 
closer to each other than the nominal resonance distance. 

Note that when the system moves away from the semi-major 
axes nominal resonant ratio, both elliptic ACR solutions tend to 
zero eccentricities (Figs. [3] and [4}. 



Higher degrees of development 

In order to evaluate the upper bound of the eccentricities for 
which the situation that we described is realistic we compare 
the results obtained for increasing degree of development of the 
Hamitonian. If the global structure of the phase space remains 
the same for higher degrees and the number, natures and posi- 
tions of ACR seem to converge we can consider that our model 
is realistic. 

We realized different tests in order to compare the structure 
of the phase space for different degrees. It seems that the degree 
of development of the Keplerian part is not significantly affect- 
ing the structure of the phase space as long as we keep the first 
four degrees terms. The main effect of the subsequent terms is to 
slightly shift the structures. However, the degree of development 
of the perturbative part is more important. 

Figures [5] and [S] show the ACR positions of the 2: 1 MMR as 
functions of G in the same case than before but with a resonant 
Hamiltonian developed up to degree 4-2, 4-3, and 4-4. Figure [5] 
show the ACR positions in the directions of cos cr,, and Fig. p\ 
show them in the directions of sin cr, . We only plot the degree 
4-2 in Fig|6]because at degrees 4-3, and 4-4, all ACR solutions 
have sin <x, =0. 

We see that around G = Go the structure of the phase space 
remains the same as for degree 4-1. However, the presence of 
degree two terms induces a significant change in the position of 
the libration center in the direction e2COS<X2. This shift is also 
observed at degrees 4-3 and 4-4. For G < G (l - 2 x 10 -3 ) 
a more complex structure appears at degree 4-2 with two new 
fixed points that diverge from the libration center in the direction 
of sincr,. Note that this structure disappears at degrees 4-3 and 
4-4. We thus interpret this structure as an artifact due to a too low 
degree of development. Between degree 4-3 and degree 4-4 there 
are no significant changes in the structure of the phase space 
but the degree 4-4 brings small corrections to the positions of 
ACR solutions. We thus conclude that for the small eccentricities 
considered here, the structure would no change qualitatively if 
we consider higher degrees of development. 

Moreover, the structure of the phase space at higher eccen- 
tricities is not important for the following discussion since we 
are interested in the very end of a tidal circularization process of 
a resonant system when planets have low eccentricities. 

The same analysis can be done for any first order MMR. For 
instance, Figs. [7] and[8]show the ACR positions of the 3:2 MMR 
in the same conditions (masses, etc.) as previously and with a 
Hamiltonian developed up to degree 4-1 (Pig jVj, 4-2, 4-3, and 
4-4 (PigjSJ. The structure of the phase space is very similar at 
low eccentricities. As for the 2:1 resonance, at degree 4-2, some 
complex structures appear for low values of the parameter G. 
However these structures are completely different between de- 
grees 4-2, 4-3, and 4-4 so they should also be consider as ar- 
tifacts. For G > Go(l - 1.5 x 10~ 3 ) the structure of the phase 
space is the same between degrees 4-3 and 4-4, and is very sim- 
ilar to the structure observed for the 2:1 MMR. This should not 
change much if we consider higher degrees of development and 
this structure seems to be common to all first order MMR. 



3. Dynamics of two resonant planets with 
dissipation 

In this section we consider the case of two planets in resonance 
in the presence of a dissipative force (tidal effect with the star, 
planet-disk interactions, etc.). In this case, G and T are no longer 
constants of the motion. However, if the dissipative force is suf- 
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Fig. 5.2:1 ACR positions in the directions of cos cr ; as functions of G. The conditions are the same as for Fig. |]but the Hamiltonian 
is developed up to degree 4-2 (left), 4-3 (center), 4-4 (right). At degrees 4-3 and 4-4, all ACR have sincH = sincr 2 = 0. For the 
degree 4-2, the positions of ACR in the directions of sin <r, are plotted in Fig. [6] 
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Fig. 6. 2:1 ACR positions in the directions of sin cr, as func- 
tions of G. The conditions are the same as for Fig. [3] but the 
Hamiltonian is developed up to degree 4-2. 



ficiently weak, their evolution is slow and we can apply the adi- 
abatic invariant theory (see |Henrard|1982| |Henrard & Lem aitre 
1983) 1. This means that the short term evolution of the system is 
still close to the one observed in the conservative case. On the 
long term, the dissipation affects the constants G and F and thus 
the parameter G. When G evolves, the phase space of the system 
evolves as presented in the previous section. 

Depending on the type of dissipation, G can either increase 
or decrease. Its evolution depends on how the dissipation affects 
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Fig. 7. 3:2 ACR positions as functions of G. The star mass is 
given by mo = m Q and planets masses by ;«i = ni2 — m®. The 
Hamiltonian is developed up to degree 4-1. Elliptic (stable) ACR 
are plotted using continuous lines whereas hyperbolic (unstable) 
ACR are plotted using dashed lines. At this degree of develop- 
ment, all ACR have sin cr 1 = sincr 2 = 0. 

the eccentricities and the semi-major axes of both planets. More 
precisely, the derivative of G is given by: 

dG 



dt 



= -A, 



n 

,2 \<?2 

'2 



(39) 
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where (e,/e,)ld, (dr/a)|d are given by the dissipation model and 
A,- are the renormalized (by F) actions. 

3.1. Migration 

In the case of a strong migration of the planets and with low ec- 
centricities, the dominant term in Eq. < (39] > is the third one which 
depend on the semi-major axes ratio evolution. For low eccen- 
tricities we have: 



-I<0. 



(40) 



Thus, if the planets undergo convergent migration ( (a/a)\ d > 0), 
G should decrease and if the migration is divergent, G should 
increase. As it is already known (see Henrard & Lemaitre 1983}, 
a resonant capture can occur only if the migration is convergent 
and the resonance is crossed with a decreasing parameter G. 

3.2. Tidal circularization 

3.2.1 . Evolution of the parameter G 

In the case of circularization of the orbits due to tides raised on 
the planets by the star, the angular momentum of each planet is 



almost conserved and the dissipation in eccentricities reads to 
lower order in eccentricities (see Papaloizou|201 



de, 



to/ 



(41) 



where t c j is a damping constant (see Eq. (82i). Since the an- 
gular momentum of each planet remains unaffected by the dis- 
sipation, the semi-major axes compensate the eccentricities de- 
creases. The evolution of semi-major axes then reads to lowest 
order in eccentricities: 



d«, 
~d7 



de, 

= lata — 
d df 



(42) 



In this case, the easiest way to compute how the dissipation 
affects the constant G is not to use Eq. ( [39] ) but to go back to the 
definition of G. Since G — G/F and G is conserved, we have: 



(43) 



From Eq. ( 42 1 we deduce the impact of the circularization pro- 
cess on A; (to lowest order in eccentricities): 



dA, 
~d7 



o2 A, 

''V 



(44) 
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Therefore, the evolution of T is given by: 



dr 

d7 



1 + 



q\ 2 A i 



,A 2 



< 0. 



d \ Pi fc,l fc,2 

And the evolution of G is governed by: 



dG 
~dT 



,Ai 



+ e 



''t c ,2 



> 0. 



(45) 



(46) 



It results that, in the case of tidal circularization of the or- 
bits, G slowly increases with time (dominant term of order 2 in 
eccentricities). Looking at Figs.[3]to[8] the position of the system 
on these graphs will slowly migrate from the left to the right, and 
eventually, the system will pass through the bifurcation (around 
G = Go). The only remaining ACR solution is the one that corre- 
sponds to the libration center but the separatrix of the resonance 
and the resonant area disappear. The motion around the only 
remaining fixed point is elliptic and this fixed point is quickly 
decreasing to zero eccentricities. Thus, the motion should be 
very well approximated by the secular problem (using a com- 
plete averaging instead of the partial one introduced to obtain 
the resonant normal form). Looking at Fig. |4j we see that the re- 
maining ACR depart from exact commensurability with a lower 
value of a (and higher value of P2/P1) than the nominal reso- 
nance value. This corresponds to planets farther away from each 
other than the nominal resonant distance. Therefore, as noticed 



by Papaloizou & Terquem ( 2010 1 the long-term effect of the tidal 



circularization is a repulsion of both planets. This is why Batygin 
|& MorbidelTI| ( [20T2] l and |Lithwick & Wu| ( |20"l2l l invoke this pro- 
cess to explain the excess of Kepler systems whose planets are 
just slightly further away from each other than nominal reso- 



3.2.2. Motion around the libration center 

The dissipation strongly affects the eccentricities (dominant 
terms of order one in Eq. (|4Ti) whereas the constants T (or 
G) are only slightly affected (dominant terms of order two in 
Eqs. (45 i,(46 1). Thus, we can consider that the damping of ec- 



centricities happens on a shorter time scale than the evolution of 
T and G. These quantities can thus still be considered as con- 
stants (in first approximation) when we take into account the 
damping of eccentricities in the equations of motion (adiabatic 
invariant theory). From Eq. ( |4T) , we deduce: 



dr 



with t c i 



Xj 

t 

T c,i 
tcj_ 

' F 3 ' 



(47) 



Now, suppose that the system lie in the vicinity of the li- 
bration center of a first order MMR whose position is noted x® 
(i = 1,2). Let us define the vector x as: 



x - 



x 2 

X\ 

K x 2 



(48) 



We can linearize the equations of motion (of the conservative 
case) around the fixed elliptic point x°: 



dx 
d7 



iAj.o (x - x°) . 



Equations of motion of the dissipative problem are thus sim- 
ply given (to first order in eccentricities) by: 



x ~ \A x d (x - - 6Bx, 
with (Eq. (|47j»: 



(50) 



6B = 



Thus: 



^ ) 

TTi 
0^0 

h, 



(51) 



x = (\A-5B) (x -x 1 ), 
with 

x l = (L4 - 6By l iAx° = (i + iA^SBf 



(52) 



(53) 



Since SB is a pertubation of L4, we can make the approximation 
(to first order): 



(i + LT'cJfi) ' * (i-LT^fi). 
Hence: 



x l ~ (i-tr'cjfi)* = x° 



iA^SBx . 



(54) 



(55) 



The dynamics around the libration center is thus modified 
by the dissipation in two different ways. First, the position of 
the fixed point is slightly changed by the offset A~ l SBx° in the 



imaginary directions (see Eq. (55 1). Secondly, the matrix giving 



the motion around the fixed point is modified (see Eq. (52 1). 

Noting y = x - x l , we can apply directly the method de- 
scribed in |Laskar et al. ( |2012| l to y. This method highlights 
the effects of the pertubation on the eigenvectors (i.e. the di- 
agonalizing linear transformation) and on the eigenvalues of 
the system. If we note So the diagonalizing matrix and Do = 
diag(gi,g2, -gi, -gi) the diagonalized matrix of the conserva- 
tive case such as: 



y 

Do 
u 



S u, 
So'ASo, 



(56) 
(57) 
(58) 



where u is the vector of eigenmodes and g\, gj are the eigenval- 
ues. The diagonalizing matrix and the diagonalized matrix in the 
dissipative case are given by: 



S = SoCI + WSO, 
D = iD -6D u 



(59) 
(60) 



The small change 6S 1 of the eigenvectors (see Laskar et al. 2012| l 
does not introduce a major change in the dynamics but can be 
seen as a corrective term as well as the correction on the position 
of the fixed point. However, the pertubation SD\ of the eigenval- 
ues is much more important. Noting 5D\ — diag(yi, 72, 71,72), 
we have: 



y i ^(S 1 6BSo\ i , 



(61) 



and these coefficients are real and positive. Finally the diagonal- 
ized equations of motion now read: 



(49) it = diag(ig! - 7i,ig 2 - 72, -igi - Ju-igi - yi)U. 



(62) 
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Thus, the dissipation introduces negative real parts {-yi) in 
the eigenvalues of the system. Therefore, all the eigenmodes will 
be damped on time scales given by these coefficients and the 
fixed point is no more elliptic but attractive. The main short-term 
effect of the dissipation is to induce an attraction of the system 
towards the libration center. 

To sum up, the tidal circularization process induce two main 
effects for planets that are initially resonant. On the short-term, 
planets tend to reach the stable ACR solution corresponding to 
the libration center and whose position is slightly modified by 
the dissipative terms. On the long-term, planets tend to leave the 
resonance by moving away from each other (the parameter G 
increases and the planets follow the stable ACR outside of the 
resonance). 

Note that, even outside of the resonance, the attractive ACR 
solution continues to exist and is not exactly at zero eccentric- 
ities. We will focus in the next section on the behavior of the 
resonant angles outside the resonances. 



4. Secular motion at very low eccentricities 

In this section we consider a conservative system of two planets 
around a star (on the same plane) that are far from any resonance 
and which have very low eccentricities. More precisely, we sup- 
pose that eccentricities have already been damped by the dissi- 
pation and we study the dynamics of the system after the dissipa- 
tion process. We are in the field of application of the secular nor- 
mal form. We start our study with the Poincare rectangular as- 
trocentric coordinates: (Ai, A,-, )>/, — iy ; ) for each planet (i = 1,2) 



(e.g. Laskar & Robutel 1995 1. Note that, for the simplicity of 
notations, we omit the hats even if these variables are not renor- 
malized anymore by T. y\,y2 are the usual Poincare rectangular 
coordinates that are noted x, x' in |Laskar~& Robutel ( 1995| l: 



Va7^i- Ji^f^ = V 



(63) 



The Hamiltonian can be developed in power series of y,, y t 
and in Fourier series of the mean longitudes /}, : 



<W = <H (A) + <Hi(A,A,y,y), 

Wi = J] h m MVy?y?y?y?e KkuM \ 

k,m,m 

where the D' Alembert rule reads: 



ki + £2 + "ii + m2 — mi — ni2 = 0, 



(64) 
(65) 



(66) 



and the coefficients hk, m jn are functions of Ai, A2 which can be 



expressed in terms of Laplace coefficients (see Laskar & Robutel 
[1993) . 

To first order of the planets masses, the secular normal form 
is obtained by averaging the perturbative part of the Hamiltonian 
over the mean longitudes. This is performed by introducing 
a change of coordinates close to the identity. Let us note 
(At, Ai,yi, — iy t -) the new coordinates, < H the new Hamiltonian, 
and W\ the generating Hamiltonian of the transformation. By 
definition, if we note 'Ho and 'Hi the Keplerian part and the first 
order part (in planets masses) of the new Hamiltonian, the trans- 
formation reads (to first order in planets masses): 



•Ho = 7f , 

#j = <H X +{W u ( Ho) 



(67) 
(68) 



where the Poisson brackets are noted with braces. The secular 
normal form (at the first order) is obtained by imposing: 



Wi=Wi + {Wi,flb} = CHi>. 



(69) 



Eq. ( |69] i is commonly called the homological equation and its 
solution is given by: 



Z^k,m,m 
_ i(mk\ + 



%,m,m(A) _ m 3BI, -mi-™! Hki'At+kih) 



fct(0,0),m,m 



"2^2) 



Ji y\ y 2 yi e 



(70) 



where m, n 2 are the unperturbed Keplerian mean-motions of the 
planets. 

By construction the secular Hamiltonian reads to first order 
of the masses: 



(71) 



Ai and A2 are constants of the motion (first integrals) and 
y\ — yi — is a stable fixed point around which the equations of 
motion can be linearized (Lagrange-Laplace theory). The change 
of variables is constructed to be close to the identity, but the re- 
version to the initial set of coordinates reintroduces short period 
terms. This reversion reads up to the first order in planets masses: 



Ai = Ai + iW u Ai} = Ai + 
Ai = A,+{Wi,A,} = A,- 

yi = yi + {Wiji} =5>; + i 



dX~ 

W, 

dWi 



1 



(72) 
(73) 
(74) 



In general, this reversion to the original set of variables only in- 
troduces small corrections to the dominant secular terms. But if 
secular eigenmodes are totally damped (i.e. yi = yi — 0), the 
only terms that are nonzero in the expressions of y\ , j2 are those 
corresponding to first order mean-motion commensurabilities: 

yi = i^L = V fe (-^i).(0.Q).(i,Q)( A ) c i(( P+ m,-pl,) ^ (?5) 



3'2 



.dW^ = y ^(-p,p+l),(0,0),(0,l)(A) ci((p+1)32 _ p3l ) 



Now, if we consider a system which is near a first order 
mean-motion commensurability (p + l):p but still outside it, the 
divisor (p + 1 )«2 -pnj is smaller than the others (but not consider- 
ably small) and the corresponding term dominates the evolution 
of y\ and y%- Thus y\ and y2 are of the form: 



a h (-P,P+ 1 ),(0,0),( 1 ,0) ( A ) e i((p+ l)A 2 -p\) 

n 2 (p + 1) - nip 

„ _ h (-p,P+ 1 ),(0,0),(0, 1 ) ( A ) ^i((p+ l)A 2 -pi 1 ) 

n 2 (p + l)-n l p 



(77) 
(78) 



In terms of arguments of periastron (see Eq. ( 63 1), this means 
that: 



(p+l)A 2 -pAi +e h 



(79) 



where e,- = or n, depending on the sign of the amplitudes in 



Eqs. (77 1, (78 1. Of course, in this situation the resonant angles 
°~\ — (p + 1)^2 — pA\ — to\ and cr 2 = (p + 1)A2 — pA\ — a>2 will 
librate (around or n). Figure |9] illustrates this phenomenon of 
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Fig. 9. Illustration of the artificial libration of the resonant an- 
gles. We plot the evolution of the complex variables e\e la> ' (Al, 
A2) and eie 10 " 1 (B 1, B2) in the case of a strong secular eigenmode 
(Al, Bl) and in the case of a severely damped one (A2, B2). For 
the sake of simplicity, we suppose that there is only one secu- 
lar eigenmode (of frequency g), but the reasoning would be the 
same with more. We also suppose that there is only one high fre- 
quency term, corresponding to the nearest first order resonance, 
and whose frequency is v = (p + X)ni-pn\ . The secular contribu- 
tion is plotted in green in all graphs whereas the high frequency 
term contribution is plotted in blue. The red curves correspond 
to the sums of both contributions. When there is a strong secular 
eigenmode, the high frequency term only brings small perturba- 
tions around the secular evolution of the argument of periastron 
a>\ (Al) and the resonant angle cr\ (Bl) which circulates. When 
the eigenmode is damped, the evolution of a>\ (A2) and <j\ (B2) 
are dominated by the high frequency v. The amplitude corre- 
sponding to v does not change between Al, Bl and A2, B2 (the 
scale changes) but in B2 o~\ appears to librate because the circu- 
lation center (marked with a blue dot in Bl, B2) is not at zero. 



artificial libration. It is important to notice that contrary to the 
resonant case, (p + 1)A.2 - pA.\ and are dominated by short 
periods (high frequencies). In the resonant case, we always have 
(p + 1)«2 — pn\ ~ 0, thus (p + 1)^2 - pA\ has a long period and 
a>i are dominated by secular eigenmodes, which also have long 
periods. This is the reason why the classical secular averaging is 
valid in this case. 

Thus, in the case of very low eccentricities, when eigen- 
modes are almost completely damped by dissipation (Figj9] 
A2,B2), the fact that the resonant angles are librating (Figj9] B2) 
does not mean that the system is resonant. It is just a geomet- 
rical effect due to the fact that the circulation center is not at 
zero. More precisely, the amplitude of the considered argument 



(Eqs. (77 1,(78 i ) is larger than the amplitude of the proper mode. 



5. Numerical simulation 

In order to confirm the different behaviors of resonant systems 
with dissipation that we highlighted in section [3] we ran two 
numerical simulations, one in the case of the 3:2 MMR and the 
other for the 2: 1 MMR. 



5.1. Simulation SI, 3:2 MMR 

The first numerical simulation concerns the two innermost plan- 
ets of the GJ581 system which are near the 3:2 MMR (e.g. 
Papaloizou 201 1). Our simulation is comparable to the one pre- 
sented in |Papaloizou| ( |201 1) . The star mass is mo = 0.3 lm , 
and planets masses are set to m\ - 1.94m e and mi - 15.64m e . 
The initial semi-major axes are set to fli = 0.11 AU and 
«2 = 0.15 AU. Initial eccentricities are set to e\ = 0.01, and 
e2 = 0.001. At the beginning of the simulation both planets un- 
dergo a migration process due to interactions with a disk. We 



adopted the same migration time scale as |Papaloizou (2011 



f mig ,,= 4.375 xl0 5 ^yr. 

fill 



(80) 



Because the outer planet is more massive, its migration is more 
efficient and the system undergo a convergent migration. During 
this migration phase the time scale of the circularization is set 
to be 100 times shorter than the migration time scale (see Lee & 
|Peale|2002| i: 



t c ,i = 



l mtgj 

Too ' 



(81) 



At t — 10 kyr, the migration stops (the disk is gone) and the 
only dissipative force that remains is the tidal interaction with 
the star. The timescale of this circularization process is given by 
(still following [Papaloizou| ( |201 l| l): 



where Q' is a parameter of the tidal dissipation model. Q' is set 
to 1 .5 in this simulation in order to reduce the computation time 
required to see the effect of the circularization (see Papal oizou] 
|20 lip . A realistic value would be several order of magnitude 
greater. This means that the time scale of the evolution of the 
system during the circularization process in our simulation is 
much shorter than the realistic one. We perform a direct n-body 
integration of the system using the ODEX integrator (see |Hairer| 



(2/3) 



v 6.5 



by: 



F d =~ 



et al. 20 1 ), with the dissipative force acting on each planet given 



1 dr, 



(83) 



Figure 10 shows the evolution of semi-major axes, and ec- 
centricities of both planets as well as the evolution of the param- 
eter G at the beginning of the simulation (up to 20 kyr). During 
the first 10 kyr of the simulation, the planets undergo convergent 
migration and G decreases (Fig. [10] bottom in gray). The sys- 
tem enters very quickly in the 3:2 resonance and we can see on 
Fig. 10 middle, that the semi-major axes stay locked in the reso- 
nant ratio. As explained by |Lee & Peale| ( [2002| l, the eccentricities 
are excited by the resonance until they reach an equilibrium state 
that depends on the strength of the circularization process which 
comes with the migration (Fig.[l0] top). 

At t — 10 kyr, the migration stops and the only remaining 
dissipative force is the tidal circularization of the orbits. The long 
term evolution of G is shown in Fig. [TT] As expected, we see on 



Fig. 1 1 that the effect of the tidal circularization is an increase of 



G that eventually exceed Go. The increase is slower and slower, 
and at the end of the simulation G is almost constant. 

We can follow the evolution of the system on graphs similar 
to those of section [2] Figure 12 gives the positions of 3:2 ACR 
solutions for the system considered in the simulation (masses), 
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Fig. 10. Evolution of eccentricities (top), semi-major axes (mid- 
dle), and the parameter G (bottom) during the first 20 kyr of 
simulation S 1 . The red curves correspond to planet 1 . The green 
curves correspond to planet 2. The gray part of the G curve (bot- 
tom) corresponds to the migration phase whereas the black part 
corresponds to the tidal circularization phase. 
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Fig. 11. Long term evolution of the parameter G during simula- 
tion SI. 



superimposed over the successive positions of the system on 
these graphs. The colors of the fixed points are the same as in 
section [2] In particular, the center of libration of the resonant 
area is plotted in red. The gray dots corresponds to positions of 
the system during the first 10 kyr (migration phase), whereas 
the black dots corresponds to the circularization phase. Since G 
decreases during the migration phase (in gray) the system goes 
from the right to the left of the graphs. On the contrary, during 
the tidal circularization phase (in black) G increases and the sys- 
tem goes from the left to the right of the graphs. 
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Fig. 12. Superposition of the successive positions of the system 
in simulation S 1 over the 3:2 ACR positions as functions of G. 
ACR positions are computed using our model with the parame- 
ters of the considered system (masses). The colors of ACR are 
consistent with graphs of section [2] but we inverted continuous 
and dashed lines in order to improve the visibility of the simula- 
tion dots. The gray dots correspond to the migration phase and 
go from the right to the left of the graph whereas the black ones 
correspond to the tidal circularization phase and go from the left 
to the right. 



At the very beginning of the simulation, the system is not 
yet in the resonance. On Fig. [12] the system lie at the right of the 
bifurcation and there is only one fixed point. The phase space is 
similar to the one of Figs. [2} E2,F2. Then, due to the convergent 
migration (gray dots) the system passes through the bifurcation 
and the phase space looks like the one of Figs. [2] A2 to D2. 
The system undergo oscillations around the libration center (red 
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Fig. 13. Superposition of the successive positions of the system 
in simulation S 1 over the 3:2 ACR positions in the plane (G, a). 
ACR positions are computed using our model with the parame- 
ters of the considered system (masses). The colors of ACR are 
consistent with graphs of section [2] but we inverted continuous 
and dashed lines in order to improve the visibility of the simula- 
tion dots. The gray dots correspond to the migration phase and 
go from the right to the left of the graph whereas the black ones 
correspond to the tidal circularization phase and go from the left 
to the right. 



dot). Finally, due to the circularization process (black dots), the 
system passes a second time through the bifurcation. Thus, the 
system leaves the resonance and the phase space is again similar 
toFigs.[2)E2,F2. 



We see in Fig. 12 that the system follows very well the ACR 
solution corresponding to the libration center both during the mi- 
gration and the circularization phases. We also see that the tidal 
dissipation induces a decrease of the amplitude of oscillations 
around this fixed point. This corresponds to the expected damp- 
ing of eigenmodes around the libration center (section^. 

Figure 13 shows the evolution of the semi-major axes ratio 
(and periods ratio) as a function of G for the simulation and for 
the fixed points of our model. We see that during the migration 
process (in gray), the semi-major axes ratio (a) increases until it 
reaches the resonant ratio. During this migration phase, a does 
not oscillate much, but when the system enter the resonance, a 
begins to oscillate around the libration center position. During 
the circularization process (in black), the system follows the li- 
bration center position, the oscillations are damped and the semi- 
major axes ratio decreases (the periods ratio increases). Thus the 
system quits the resonance with a periods ratio larger than the 
resonant one, as expected. 

Note that when G increases, the eccentricities quickly tend 



to zero (Fig. 12 1 and the circularization process is less and less 
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effective. This is why G increases slower and slower (Fig 
This means that the system cannot depart very far away from the 
MMR on a finite time (see appendix [B] for an estimation of the 
asymptotic evolution of G). This is even more true for a real sys- 
tem, since in this simulation the circularization process is several 
order of magnitude more efficient than in reality. This mecha- 
nism explains why Kepler systems lie just slightly farther away 
from the MMR. 

Figure [14] shows the evolution of the resonant angle <t\ (in 
red) and the difference of periastrons Aco = <x 2 - lt\ (in green). 
We see that both resonant angles and Aco librate during the whole 
simulation even if the system is outside the resonance at the end 
of the simulation (see Figs. 12 and 13 1. <T\ librates around 0, 
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Fig. 14. Evolution of the resonant angle <t\ (red) and the differ- 
ence of periastrons Aw = 1x2-0-! (green) in simulation SI. 




Fig. 15. Evolution of the eccentricities (top) the semi-major axes 
(middle) and the parameter G (bottom) during the first 20 kyr 
of simulation S2. The red curves correspond to planet 1. The 
green curves correspond to planet 2. The gray part of the G curve 
(bottom) corresponds to the migration phase whereas the black 
part corresponds to the tidal circularization phase. 



5.2. Simulation S2, 2:1 MMR 

We also ran a simulation of the same system but in the case of the 
2: 1 MMR. The masses of the star and the planets are the same 
but initial semi-major axes are set to ci\ = 0.11 AU, and «2 = 
0.18 AU. Initial eccentricities are set to e\ — 0.002, and e 2 - 0. 
All other parameters are the same as for the first simulation. 



while <r 2 and Aco librate around n. This could have also been 
deduced from Fig. [12] 



Figure 15 shows the evolution of eccentricities, semi-major 
axes and G during the first 20 kyr of the simulation. Their be- 
haviors are very similar to those of the simulation S 1 (Fig.[T0|. 
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Fig. 16. Long term evolution of the parameter G during simula- 
tion 52. 



Note that as for the resonance 3:2, G decreases during the migra 
tion phase (in gray). We plot in Fig 16 the long term evolution of 
G which increases due to the circularization process. Figures 17 
and 18 show the superposition of the fixed points of our model 



with the successive positions of the system. Figure 17 gives the 
eccentricities as functions of G and Fig. 18 gives the semi-major 
axes ratio (as well as the periods ratio) as a function of G. As for 
the 3:2 resonance, in the case of the 2:1 resonance, the system 
begins outside the resonance (on the right of the graphs) and the 
phase space corresponds to the situation of Figs. [2] E1,F1. Due 
to the convergent migration of planets the system goes from the 
right to the left (gray dots) and passes through the bifurcation 
following the red ACR solution. At this point the phase space 
is similar to those of Figsj2] Al to Dl. During the circulariza- 
tion phase, the system goes from the left to the right in Figs.[T7] 
18 (black dots). The system passes again at the bifurcation and 
leaves the resonance (phase space of Figs. [2j E1,F1). Shortly 
after the system quited the resonance, the amplitude of oscilla- 
tions around the fixed point suddendly increases in the direction 
of e 2 . This is due to the crossing of a 2:1 resonance between the 
two proper modes of the motion around the fixed point (see ap- 
pendix|C]for more details). However the dissipation damps again 
the oscillations after this event and the system keeps following 
the ACR solution. 

We observe on Fig.[l9]that the first resonant angle cr 1 keeps 
librating around during the whole simulation whereas Aoj 
and o"2 circulate almost all the time. More precisely, during the 
first 0.5 Myr of the simulation, Aoj and cr 2 librate around 0. 
Then, these angles circulate during a short laps of time around 
t = 0.5 Myr. After this event, they begin to librate around ji. 
Actually, the system is still in resonance and oscillates around 
the libration center but the position of this fixed point crosses 
the value e 2 = for G/G - 1 = 5 x 10~ 4 (see Fig. [IT]) around 
t = 0.5 Myr. The dynamics does not change and this event is 
only geometrical. Thus, even if the system is clearly in the res- 
onance area, very close the the libration center, we observe a 
circulation of cr 2 . Around t = 1.4 Myr, Aw and cr 2 begin to 
circulate (Fig. 19 1. This corresponds to the crossing of the res- 
onance (see appendix [C| that we observe in Fig.[l7] This event 
increases the amplitude of oscillations around the fixed point and 
due to this higher amplitude the system passes through <? 2 = at 
each oscillation. This is why we observe a circulation of Aa> and 

0"2- 

Note that both simulations show that it is possible to observe 
an oscillation of the resonant angles outside of resonances (a~i 
and o"2 for S 1, and only cr\ for S2). On the opposite, in simula- 
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Fig. 17. Superposition of the successive positions of the system 
in simulation S2 over the 2:1 ACR positions as functions of G. 
ACR positions are computed using our model with the parame- 
ters of the considered system (masses). The colors of ACR are 
consistent with graphs of section [2] but we inverted continuous 
and dashed lines in order to improve the visibility of the simula- 
tion dots. The gray dots correspond to the migration phase and 
go from the right to the left of the graph whereas the black ones 
correspond to the tidal circularization phase and go from the left 
to the right. 



tion 5 2 we observe a circulation of <x 2 when the libration center 
position crosses the axis <? 2 = while the system is clearly inside 
the resonance and oscillates very close to the libration center. 
Therefore, it appears that resonant angles have to be considered 
with great caution and cannot always be used to distinguish be- 
tween truly dynamical effects and simple geometrical effects. 
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Fig. 18. Superposition of the successive positions of the system 
in simulation 52 over the 2: 1 ACR positions in the plane (G, a). 
ACR positions are computed using our model with the parame- 
ters of the considered system (masses). The colors of ACR are 
consistent with graphs of section [2] but we inverted continuous 
and dashed lines in order to improve the visibility of the simula- 
tion dots. The gray dots correspond to the migration phase and 
go from the right to the left of the graph whereas the black ones 
correspond to the tidal circularization phase and go from the left 
to the right. 
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Fig. 19. Evolution of the resonant angle ct\ (red) and the differ- 
ence of periastrons Aa> = o"2 - en (green) during the first 10 Myr 
of simulation S2. 

6. Conclusion 

In this work, we presented a study of planar resonant planetary 
systems in the conservative case and in presence of a dissipative 
force. Our main interest was to understand the dynamics at the 
end of a circularization process in resonant systems. We used a 
completely analytical model developed in power series of eccen- 
tricities which is well suited for the study of these low eccentric- 
ity systems. Before introducing the dissipative force in the model 
we characterized the dynamics in the conservative case. In par- 
ticular, we highlighted the fact that apsidal corotation resonances 
(ACR) are a powerful tool to understand the global dynamics of 
a system. Then, we showed that the introduction of a dissipative 
f orcein a resonant system has two main effects. On the short- 
term, the system is attracted toward the libration center if it ini- 
tially relies in its vicinity. On the long-term, the system tends to 
follow this ACR solution outside the resonance and both planets 
tend to move away from each other. These two mechanisms are 
well illustrated and confirmed by the results of two simulations 
(one for the 2:1 MMR and the other for the 3:2 MMR) of the 



innermost planets of the GJ581 (see section [5J. Since the ACR 
solution do not correspond to zero eccentricities even far from 
the resonance, it is possible to have resonant angles to oscillate 
outside of resonances. However, we showed that in this case the 
motion is completely characterized by the secular problem and 
that the fact that resonant angles appear to librate only means that 
the secular eigenmodes are (almost) totally damped. The impor- 
tant fact is that the nature of the motion is the same as when the 
eigenmodes are not damped and the separatrix of the resonance 
does not exist anymore in this region of the phase space. Thus, 
it is inappropriate to consider such systems as resonant ones. 
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Appendix A: Computation of the Hamiltonian 
coefficients 

In this section we explain in more details the computations of 
coefficients K, S , and R of the resonant normal form (Eq. pi} ) 
up to any order. 

The Keplerian part is given by: 



1 We consider here tidal effects, but it is clear that other mechanisms 
could result as well in similar dissipative effects. 



2A 2 



(A.l) 
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with: 



ficients can be expressed as functions of a and Laplace coeffi- 
cients (fc'/V)): 



q 

A 2 = l-[l + i|Ai = (l + £)(G+D)--. 



(A.2) 
(A.3) 



K coefficients are thus simply given by the Taylor series of 
1/(1 + x) 2 where x = T>j{G - 1) for the inner planet and 
x - T>l (G - p/(p + q)) for the outer one. 

For the pertubation, the first step is to determine which in- 
equalities k\A\ + k 2 A 2 have to be computed. The secular part is 
given by the inequality k\ - k 2 - 0. The resonant terms are given 
by inequalities k\ = -kp, k 2 = k(p + q), with -d/q < k < d/q 
where d is the chosen degree of development. Each inequality 
is computed as a power series of x\, X2 by using the algorithm 
presented in Laskar & Robutel ( 1995 ). For instance, in the case 
of the 2:1 resonance, the perturbative part of the Hamiltonian 
reads: 
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where 1ri\ >( j is the direct part of the development, and , 7Yi ) j is the 
indirect part. The direct part is given by: 



+ rf\x 1 +X l ) + / 2 ]l(X 2 +X 2 ) 

+ rfiXl+t^fiXl+xl) 
+ tto(XiX 2 + X X X 2 \ 



(A.5) 



r (2) 
'12 



_L3 



> 69 
-a + — a + 



36 



40 

29 59 



'3/2V 



5^)6® (or) 



24 



+ ir 2 + 2o + 2o ffZ + y^T^) 

As for the Keplerian part, we substitute Ai and A 2 by their 
power series in T> each time they appear in the expression of 
"Hi . For the Laplace coefficients we first need to develop them in 
power series around the nominal resonant semi-major axes ratio 
which is defined by: 



and the indirect part reads: 



1/3 



p + q 



2/3 



(A.8) 



'Hii = r 2 1 }(X 2 + X 2 ), 



(A.6) 



Then the deviation 6a = a-a r to the nominal value is substituted 
by its power series in D. Of course the degrees of development 
of all these series are adjusted in order to be consistent with the 
desired degree in eccentricities (d). 



where we note X, 



Laskar & Robutel 



1995 ). The coefficients S , R of Eq. pT} correspond to coeffi 
cients s, r appearing in Eqs. (A.5i, ( |A.6[ ) but with a renormal- 
ization due to the use of x\, x 2 instead of X\, X 2 and the factors 
in front of the direct and indirect parts in Eq. dA4|. s, r coef- 



Appendix B: Asymptotic evolution of G 

In this section we show how to estimate the asymptotic evolution 
of G. This allows to determine the position of the system in our 
graphs as a function of time. We have seen that the evolution 
of G during the tidal circularization phase is slower and slower. 
Moreover, the escape from the resonance is very quick and the 
evolution of the system after this escape is very slow. Thus, the 
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initial position of the system in the resonance does not influence 
much its final position (outside the resonance) after a long time. 

The evolution of G during the circularization phase is gov- 
erned by Eq. ( po*) . Since G remains very close to Go, and A, stay 
very close to A,- >r , we can substitute these values in the expres- 
sion of G. However, e, cannot be considered as constants since 
the ACR solution tends to zero eccentricities when G increases. 
Thus we have to compute e-, as functions of 6G = G-Gq with the 
help of Eqs. (37 1, (38 i. In the asymptotic evolution, eccentrici- 
ties are very low and the ACR position can be well approximated 
with lowest order terms : 



* K™x x +Rf, 
* K™X2+S$\ 

Thus the ACR solution is simply given by: 



R 



x\ = - 



KO-r 



R 



(B.l) 
(B.2) 

(B.3) 
(B.4) 

(B.5) 



K (2 \ Rf\ and K£> are functions of G but R™, and R%> are almost 
constants and can be approximated with their values for G = Gq 
(see appendix [A): 



and <?, are given by: 
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On the contrary, in first approximation is proportional to SG: 

\2 ,,2o3" 
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Finally, this means that: 
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The asymptotic evolution of G is thus given by: 

5G oc t x l\ 

The eccentricities evolve as: 



e; OC t 



■1/3 



(B.8) 



(B.9) 



(B.10) 
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(B.12) 



The semi-major axes ratio asymptotic evolution is governed by: 
6aocSGoct l/i . (B.13) 



We thus find the same asymptotic law (f 1 ^ 3 ) as Papaloizou 
pOTT) , and |Lithwick & WulpOTI] ). 



Appendix C: 2:1 resonance between proper modes 
around the ACR solution in S2 



0.15 
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Fig. C.l. Superposition of the trajectory of the system 52 during 
the event (around t — 1 .4 Myr) with the energy levels curves of 
the 2:1 resonance between both proper modes around the ACR 
solution. The successive positions of the system just before the 
event are marked with red dots while the positions of the system 
just after the event are marked with green dots. 



In this section we present a brief analysis of the event hap- 
pening in S2 around t — 1.4 Myr (Fig. 



17 i. We interpret the 



increase of the amplitude of oscillations in the direction of e 2 
as a consequence of the crossing of the 2:1 resonance between 
both proper modes around the ACR solution. The linear equa- 



tions of motion around the fixed point (see Eq. (58 1) correspond 
to degree two terms in the diagonalized Hamiltonian: 



diag,2 - ~g\U\U\ - g 2 M 2 M 2 



(CI) 



The resonance appears in this Hamiltonian at degree three with 
the term of the form: 



(C2) 



where p us a constant depending on the parameters and initial 
conditions. Let us note = Vt^e lt,f . We can average over all 
other angles than the resonant one (v\ — 2i/ 2 ) with a similar pro- 
cedure than the one used in section|2] The system is then reduced 
to one degree of freedom and we have a new constant of the (av- 
eraged) motion: 



U = 2Ui + U 2 . 



(C3) 



The energy levels curves for the values of the constants G and 
U taken by the system of simulation 52 at the moment of the 
event (around t = 1.4 Myr) are shown in Fig. |C.1| superimposed 
with the trajectory of the system. We see that before the event, 
the proper mode 2 is almost totally damped compared to proper 
mode 1 (red dots in Fig. |C.1[ ), whereas just after the event, this 
proper mode gained some amplitude (green dots in Fig. C. 1 1 and 



the system follows very well the expected motion. The resonant 
angle v\ - 2i/ 2 does not enter in libration in the simuation S 2 but 
this may be due to the very strong dissipation present in this sim- 
ulation and which induces an evolution of the constant G quicker 
than expected for a real system. Thus the frequencies g\, g 2 and 
the phase space evolve very quickly and the system crosses the 
resonance before it has time to enter in the libration area. The 
important fact is that due to the crossing of the resonance, the 
proper mode 2 increases its amplitude while the proper mode 1 
amplitude's decreases (U = 2U\ + t/ 2 stays constant). However, 
since the proper mode 1 has initially a much higher amplitude, 



17 



J.-B. Delisle et al.: Dissipation in planar resonant planetary systems 



this decrease is imperceptible. Moreover, since the diagonaliz- 
ing matrix is dominated by diagonal terms, the evolution of e\ 
is dominated by the proper mode 1 and the evolution of e-i by 
the proper mode 2. Thus, e\ is only weakly affected by the event 
while e2 is strongly affected. 
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